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ABSTRACT 


/ 


-We  disctwr finite  difference  methods  for  partial  differential  equations  on  polar  and 
spherical  coordinate  systems.  The  distinctive  feature  of  these  coordinate  Systems  is  the 
coordinate  system  singularity  at  the  origin.  We  show  how  to  accurately  and  conveniently 
determine  the  solution  at  the  origin  for  both  scalar  and  ‘or  fields.  \8?e  also  discuss 
the  Fourier  method  to  approximate  derivatives  with  respect  to  the  angular  variable  in 
polar  coordinates.  Computational  examples  are  presented  illustrating  the  accuracy  and 
efficiency  of  the  method  for  hyperbolic  and  elliptic  equations,  and  also  for  the  computation 


of  vector  fields  at  the  origin. 
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SIGNIFICANCE  AND  EXPLANATION 


Many  finite  difference  computations  are  done  for  regions  with  polar  or  spherical 
coordinate  systems.  The  determination  of  the  solution  variables  at  the  origin  of  such  coor¬ 
dinate  systems  has  been  a  source  of  much  confusion.  This  is  especially  true  for  calculations 
of  vector  fields  uefined  on  polar  or  spherical  grids.  In  this  paper  we  show  how  to  accurately 
and  easily  calculate  the  variables  at  the  origin  of  such  coordinate  systems. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive  summary  lies 
with  MRC,  and  not  with  the  authors  of  this  report. 


FiNlTE  DIFFERENCE  METHODS  FOR  POLAR  COORDINATE  SYSTEMS 

John  C.  Strikwerda  and  Yvonne  Nagel* 

1.  Introduction 

In  this  paper  we  consider  the  use  of  polar  and  spherical  coordinates  with  finite  dif¬ 
ference  methods,  determining  how  to  achieve  accurate  results  with  convenience.  Although 
the  use  of  finite  difference  methods  with  polar  coordinates  is  not  at  all  new  there  are  sev¬ 
eral  features  of  their  use  that  are  not  well  known  among  numerical  analysts,  computational 
scientists,  and  engineers.  In  particular,  the  accurate  treatment  of  vector  fields  with  polar 
coordinates  is  not  widely  known. 

It  is  the  aim  of  this  paper  to  bring  together  the  pertinent  information  and  present 
it  in  an  organized  way.  As  such,  this  paper  presents  few  new  ideas,  but  it  is  hoped  that  it 
will  be  a  useful  addition  to  the  literature  on  numerical  methods. 

Much  of  what  is  presented  here  also  applies  to  axially  symmetric  problems  in  polar 
or  spherical  coordinates;  the  common  feature  of  these  problems  is  the  singular  nature  of 
the  coordinate  system  at  the  origin. 

2.  The  Center  Formulas 

Consider  the  plane  with  a  polar  coordinate  system.  Each  point  is  determined  by 
its  polar  coordinates  {r,4>)  which,  for  points  other  than  the  origin,  is  unique  up  to  integer 
multiples  of  27r  in  <t>.  However  the  origin  has  the  coordinates  (0,<^)  for  all  angles  <f>,  and 
it  is  this  lack  of  uniqueness  in  the  coordinates  of  the  origin  that  introduces  difficulties  for 
numerical  methods.  These  difficulties  are  displayed  in  the  Jacobian  of  the  coordinate  map 
which  takes  ordered  pairs  in  jO,  oo)x  IR  to  points  in  the  plane.  The  Jacobian  vanishes 
on  {  0  }  X,  IR.  However  it  is  important  to  realize  that  this  singularity  is  present  in  the 
coordinate  map  and  polar  representations  of  functions  and  need  not  be  present  in  the 
functions  themselves.  In  this  paper  we  shall  only  consider  functions  which  are  smooth  in 
the  domain  being  considered. 

The  singular  behavior  of  the  polar  coordinate  system  at  the  origin  usually  precludes 
the  direct  use  of  finite  difference  approximations  to  differential  equations  at  that  point.  We 
consider  therefore  the  use  of  interpolation  formulas  to  accurately  determine  the  solution  at 
the  origin.  We  begin  by  considering  a  function  defined  in  the  plane,  without  considering  a 
coordinate  system. 

'Sponsored  by  the  United  States  Army  under  Contract  No.  DA  A  029-80-0-0041.  The  first  author  was  also 
supported  in  part  by  NSP  Grant  MCS-8,30fia80  and  ONR  Grant  N00014-84-K-0454. 


Consider  a  smooth  function  u  defined  in  a  neighborhood  of  a  point  P  in  the  plane. 
We  wish  to  express  u[P)  in  terms  of  averages,  u(P,p),  on  circles  of  radius  p  centered  at  P. 
We  begin  by  expanding  u  in  a  Taylor  series  in  cartesian  coordinates  with  the  origin  at  P, 


E 


t,J=0 


kUldx^dyt 


(0, 0)  +  Rjif- 


where  R^/  is  the  remainder  term.  Then, 


(2.1) 


1  f*' 

u(P,p)  =  —  I  u{pcos^,psin^)d^ 

27r  Jo 


(2.2) 


and  using  (2.1) 


u(P,p)  =  ZcyV^‘u{P)  +  Rn 


1=0 


(2.3) 


where  cj  =  1/4*(Z!)®.  Formula  (2.3),  which  is  independent  of  a  coordinate  system,  is  the 
basis  for  determining  a  function  value  at  the  origin  given  values  of  the  function  at  points 
nearby  and  given  the  differential  equation  satisfied  by  u. 


Consider  now  a  uniform  finite  difference  grid  with  grid  points  (r,,^y)  for  integers  i 
and  j  with  t  >  0  and  0  <  j  <  J  -  1  where  r,  =  lAr,  <j>j  =  jA4>  for  Ar  >  0  and  A^  =  2ir/J. 

For  a  function  u  defined  in  a  neighborhood  of  the  origin  P,  we  have 

u(/Ar)  =  «(P,/Ar)  =  7  +  0(An  (2.4) 

''  j=o 

where  ui  j  =  u{rt,<f)j)  and  m  is  a  positive  integer  whose  value  will  be  considered  in  section 
5.  Using  (2.3)  we  then  have  the  relations 

u(P)  =  u(Ar)  +  0(Ar*)  +  0(A<^”*)  (2.5) 

and 

u(P)  =  ^(4u(Ar)-u(2Ar))  +  0(Ar')  +  0(A<^"*)  (2.6) 

which  can  be  used  with  finite  difference  methods  to  determine  values  at  the  origin.  Higher 
order  formulas  can  be  obtained  by  similar  means. 

3.  The  Laplacian  with  Polar  Grids 


When  the  differential  equation  being  solved  involves  the  laplacian  operator  then 
formula  (2.3)  can  be  used  to  a  special  advantage.  We  consider  as  examples  the  Poisson 
equation 

V*u  =  /  (3.1) 

and  the  wave  equation 


=  V*u 


on  a  disk  of  unit  radius. 


For  the  Poisson  equation  (3.1)  consider  the  semi-discrete  finite  difference,  approxi¬ 


mation 


r.  Ar 


Ar 


Ui{<t>)  -  u,- 
Ar 


where  we  discretize  only  the  radial  direction.  Employing  (2.3)  we  have  at  the  origin 

Ar* 

uo  =  u(Ar)  -  ^V*u(0)  +  0(Ar<) 

4 

or 

Uo  =  u(Ar) - — /(O) -I- O(Ar^).  (3.4) 

4 

This  formula  maintains  the  second-order  accuracy  of  the  scheme  and  is  easy  to  use.  How¬ 
ever,  even  whep  the  equation  being  solved  involves  the  laplacian,  (2.6)  may  be  more  ac¬ 
curate  or  convenient  to  use  than  formulas  such  as  (3.4).  Formula  (3.4)  has  been  used  by 
Swarztrauber  and  Sweet  jl973]  for  solving  the  Poisson  equation  in  a  disk,  and  by  Swarz- 
trauber  [1974]  for  the  Poisson  equation  on  a  sphere. 

For  the  wave  equation  (3.2)  we  also  consider  a  semi-discrete  approximation  in  which 
only  time  and  the  radial  direction  are  discretized  with  the  angular  variation  continuous. 
Let  u”{<l>)  be  the  approximation  to  u(nA<,r, ,<;>).  At  the  origin  we  have 


u”  =  uJ(Ar)  -  ^Ar^V^uS  +  0{Ar^) 
4 

=  uS(Ar)-  iAr*|^[+0(Ar^) 

A  2 


+  0(Ar*)  +  (3(Ar'A(’), 

using  a  central  difference  approximation  in  time.  This  gives  the  formula  at  the  origin  as 


«;♦'  =  K  -  ‘  *  -I  (sSlA"-)  -  “o") 


which  maintains  the  second-order  accuracy  of  the  scheme.  Example  1  in  section  7  shows 
that  this  formula  gives  accurate  results.  Similar  methods  can  also  be  used  with  parabolic 
equations. 


4.  Vector  Fields  with  Polar  Grids 


In  addition  to  the  coordinate  singularity  at  the  origin,  the  polar  coordinate  rep¬ 
resentation  of  vector  fields  introduces  an  additional  difficulty.  Let  F  be  a  vector  field 
defined  on  a  domain  on  which  there  is  a  polar  coordinate  system.  The  polar  coordinate 
representation  assigns  to  each  vector  F(P)  the  component  in  the  radial  direction  and  the 
component  in  the  direction  of  increasing  angle.  This  representation  is  unique  at  all  points 
other  than  the  origin. 

At  the  origin  the  vector  F(0)  has  a  different  representation  for  each  choice  of  the 
radial  direction.  This  is  best  illustrated  using  the  mapping  between  the  polar  and  cartesian 
representations.  Let  {U,  V)  be  the  usual  cartesian  representation  of  the  vector  field  F  which 
is  uniquely  determined,  then  the  polar  representation  (u,v)  is  given  by 

u  =  U  cos  ^  -f  K  sin  ^  (4.1) 

V  =  -  f/  sin  +  F  cos 

Since  at  the  origin  the  pair  iU,V)  is  single  valued,  (4.1)  shows  the  multivalued  nature  of 
the  polar  representation. 

Using  a  polar  grid  the  vector  field  F  will  be  represented,  and  approximated,  by 
vectors  (u,/,  v,y)  at  each  grid  point  (r,,^y).  At  the  origin,  there  is  a  representation  (uo>,Vo>) 
for  each  coordinate  direction  (0,  <f>j).  For  consistency  these  representations  must  be  related 
by  the  formulas  (4.1).  That  is,  there  are  values  (tfo,Vb)  such  that 

=  Uo  cos  +  Fo  sin  (4.2) 

vo;  =  -Uo  sin<^,  +  Fo  cos  <f>j. 

The  values  of  Uq  and  Fo  can  be  obtained  by  formulas  such  as  (2.6).  For  example,  on  a 
uniform  grid  define 

1 

U(iAr)  =  ~  cos 4>j  -  Vijsm<^>j  (4.3) 

''  ;=o 

1 

F  (j  Ar)  =  7  H  sin  0,  +  Vij  cos  <(>,. 

''  j=o 

Then  Uq  and  Fq  can  be  approximated  by 

Uo  =  ^(4U(Ar)  -  t/(2Ar)) 

Fo  =  ^{4F(Ar)  -  F(2Ar)). 


(4.4) 


These  values  can  then  be  used  in  (4.2)  to  give  the  values  of  (uq,,  vq,).  Example  3  in 
section  7  demonstrates  the  accuracy  of  this  method  2is  applied  to  the  Stokes  equations. 
This  method  has  been  used  in  Strikwerda  [1984a]  and  Nagel  and  Strikwerda  [1985]  with 
excellent  results. 

For  finite  difference  grids  which  are  not  uniform  in  the  angular  variable  formulas 
(4.3)  should  be  replaced  by 


U {iAr)  =  -  Y)  u,y(sin  -  sin ^y_i)  +  v,>(cos <i>j+i  -  cos (^>-i) 

Vi=o 

1 

V (iAr)  =  ■  I  II  «iy(cos  -  cos  ^^-i)  +  t;,y(8in  -  sin  ^>-i) 


where 


0  =  2^  sin(^,+,  -  ^,). 

3=0 


The  formulas  (4.5)  and  (4.6)  are  exact  for  the  case  when  the  vector  field  has  constant 
cartesian  components  in  a  neighborhood  of  the  origin. 

5.  The  Fourier  Method 

We  now  consider  the  Fourier  method  for  the  approximation  of  derivatives  with 
respect  to  <(>.  Consider  a  periodic  discrete  function  fj  defined  on  grid  points  =  J^4> 
with  A<t>  =  2n/J.  The  object  of  both  the  finite  difference  and  Fourier  methods  is  to  obtain 
approximations  to  d//d^  at  the  grid  points.  The  Fourier  method  begins  with  the  finite 
Fourier  series  representation  of  /,,  i.e.  for  the  case  when  J  is  an  even  integer 

JI2-1 

fj  =  5o  +  (®*  sinAr^j  4-  6*  cos  +  bj  cos{~<l>j).  (5.1) 

*=i  ’ 

Note  that  cos(^^;)  =  (-1)^.  Replacing  <f>3  in  (5.1)  by  a  continuous  variable  <t>  we  can 
approximate  df fd<t>  0,1  as 


—  H  a*A:  cos  k<i>j  —  b^k  sin  k<i>j 


and  similarly 


S  sin  k<t>,  -  bkk^  cos  kcpj)  -  (-1)^^)*  bi. 

The  coefficients  a*  and  6*  are  easily  obtained  by 


2 

6*  =  -  cos  k<l>j  (5.3) 

}=0 

2  •'-» 

at  =  sin  k<j>j 

i=o 

for  0  <  A:  <  J/2,  and 

io  =  \Ef,  (5.4) 

i=o 

=  7  E(- >)’/>• 

i=o 


The  Fourier  method  has  the  advantage  that  it  gives  far  higher  accuracy  for  a  given 
number  of  grid  points  than  do  finite  difference  methods  (Gottlieb  and  Orzag  [1977]).  Al¬ 
ternatively  to  attain  a  given  accurax:y  the  Fourier  method  requires  significantly  fewer  grid 
points  than  do  hnite  difference  methods.  For  example  2  of  section  7,  finite  difference 
methods  would  require  at  least  three  times  as  many  grid  points  in  the  angular  direction 
to  obtain  comparable  accuracy. 


The  efficiency  gained  by  the  Fourier  method  over  the  finite  difference  method  for 
the  angular  variation  is  due  to  the  natural  periodicity  in  the  variable  <l> .  Spectral  methods 
can  be  used  with  the  radial  variation  but  not  necessarily  with  the  same  gain  in  efficiency, 
Gottlieb  and  Orzag  [1977]. 


Line  successive-over-relaxation  (LSOR)  can  easily  be  used  to  solve  elliptic  boundary 
value  problems  in  polar  coordinates  in  which  the  Fourier  method  is  used  to  approximate  the 
derivatives  with  respect  to  The  basic  formula  for  LSOR  as  applied  to  the  semi-discrete 
approximation  (3.3)  is 


VnAr  \  Ar  ’"5  Ar  j 


1 

^  rf  d<l>^ 


In  (5.5)  the  order  of  progression  through  the  grid  for  the  LSOR  is  in  the  direction  of 
decreasing  radius.  When  the  Fourier  method  is  used  to  approximate  the  derivatives  with 
respect  to  the  Fourier  coefficients  of  the  update,  -  u*'  can  be  easily  obtained  from 
the  coefficients  of  the  right-hand  side  of  (5.5).  That  is,  the  right-hand  side  of  (5.5)  is 
evaluated  for  each  value  of  j,  then  the  Fourier  coefficients  are  calculated.  Dividing  the 
coefficients  of  the  node  by  -(2/Ar^  +  gives  the  coefficients  of  the  update,  from 

which  the  update  is  determined  at  each  value  of  j.  This  method  is  used  in  examples  2  and 
3  of  section  7. 


6 


6.  Quadrature  Formulas 

The  approximation  of  integrals  by  sums  arises  in  several  contexts  in  the  use  of  finite 
difference  methods  on  polar  grids.  As  we  have  seen  the  approximations  at  the  origin  (2.5) 
and  (2.6)  use  integrals  in  <j>  at  various  values  of  r.  Also,  the  accurate  determination  of 
integral  quantities  over  the  domain  requires  quadrature  formulas  in  r  and  <i>. 

We  begin  by  considering  integration  in  the  angular  variable  only.  We  consider  a 
27r-periodic  function  /(^).  We  first  consider  the  error  resulting  from  approximating  the 


integral 

/*  Sir 

by  the  sum 

/  md<f> 

Jo 

(6.1) 

j-i 

}=0 

(6.2) 

where  =  2‘k/J  and  4>j  —  By  the  theory  of  Fourier  series  we  have 


m  =  E  On*’""  (6.3) 

n=-oo 

where 

(6-4) 

Thus 

i;/W))A*=  E  E  E  ‘kj. 

j=0  }~0  n=-oo  k—~oo 

J-l 

since  ^  e"'*’  vanishes  unless  n  is  a  multiple  of  J.  The  integral  (6.1)  is  precisely  co  thus 
;=o 

the  error  in  the  approximation  (6.2)  is 

00 

27r  ^  Ckj. 

k=-oo,k9«0 


By  the  definition  of  the  c„  in  (6.4)  we  have 

an=  :^(*n)-"‘  (<^)c-"^d<^ 

Jtt  Jo 

if  /  is  m  times  differentiable.  Thus 


ia„l  <  CM'”', 


for  some  constant  depending  on  /.  Thus  the  error  in  the  approximation  (6.2)  is 
bounded  by 

2C^  £  \kJ\-”'  =  (6.5) 

t=i 


We  now  consider  quadrature  for  the  unit  circle  using  uniform  spacing  in  r  and  (j>. 
We  consider 


/  f  f(r,<i>)rdrd4)  =  2irf  f{r)rd<)>  (6.6) 

Jo  Jo  Jo 

where  /(r)  may  be  approximated  to  0(A<^“)  as  in  (6.2).  Using  the  trapezoid  formula  we 
have 


2n  f  f{r)rdr  =  2;r  i(7.r.  +  7.+,r.+j)Ar  +  ^?(Ar*) 

Jo  .•  n  Z 


i=0 

/-I 


=  27r  JiTjAr  +  nfjrjAr  +  O(Ar^). 


i=l 


Hence 


f  f  /(r,0)»’drd«/>=  +  (6-7) 

•'0  i=l  >=0  ^  >=0 


7.  Computational  Results 

In  this  section  we  present  results  of  computations  using  the  formulas  discussed  in 
the  previous  sections  applied  to  three  test  problems.  The  first  test  problem  is  to  solve  the 
second-order  wave  equation.  The  two  formulas  (2.6)  and  (3.5)  for  determining  the  solution 
at  the  origin  are  compared.  The  second  test  problem  is  to  solve  an  elliptic  equation  using 
the  LSOR  method  given  in  section  4  to  solve  the  discrete  equations.  The  third  test  problem 
uses  the  Stokes  equations  to  illustrate  the  use  of  the  formulas  for  vector  fields  at  the  origin. 

The  first  test  problem  was  to  solve  the  second-order  wave  equation 


Utt  =  V*u 


(7.1) 


in  the  unit  disk  for  0  <  /  <  1.  The  exact  solution  we  used  was 


u(<,i,y)  -  cos(t  -  .6i  -  .8y). 


(7.2) 


The  equation  for  the  time  advancement  is 


-r,"  =  (7.3) 

where  the  discrete  laplacian,  Vj,  is  given  by  the  left-hand  side  of  (3.3)  and  the  derivatives 
with  respect  to  (t>  are  approximated  by  the  Fourier  method.  The  formula  for  the  first 
lime-step  is  based  on  a  Taylor  series  in  time  and  is 

<  =  »“  +  Al)".)?,  +  ,  (7.4) 

where  and  were  obtained  from  the  exact  solution. 

Both  the  interpolation  formula  (2.6)  and  the  formula  (3.5)  were  used  to  determine 
the  solution  at  the  origin.  The  interpolation  formula  was  applied  using  u(Ar)  and  u(2Ar) 
at  the  given  time  level  to  compute  u  at  the  origin  for  that  same  time  level.  The  results 
of  four  test  cases  are  displayed  in  Table  1,  where  I  and  J  are  the  number  of  radial  and 
angular  grid  points,  respectively,  and  K  is  the  number  of  time  steps.  Both  the  norm  of 
the  error  and  the  error  at  the  origin  are  shown  for  each  case.  The  two  formulas  are  seen  to 
be  comparable  in  accuracy,  but  the  interpolation  formula  is  slightly  more  accurate.  This 
was  also  observed  in  all  the  other  cases  in  which  these  two  formulas  were  compared.  Since 
formulas  (2.6)  and  (3.5)  yielded  comparable  results,  we  used  the  simpler  formula  (2.6)  for 
all  subsequent  runs. 

Table  1.  Two  Center  Methods 


GRID 

FORMULA  1 

FORMULA  2 

I  J  K 

ll««rr|| 

i 

^err 

ll«rrr|| 

Cerr 

21  16  160 

1 

.5586(-4) 

.1010(-3) 

.5685(-4) 

.1069(-3) 

41  20  220 

j 

;  .1662(-4) 

.3177(-4) 

.1623(-4) 

.3 105  (-4) 

A  list  of  cases  using  formula  (2.6)  is  given  in  Table  2.  The  data  show  that  errors 
are  relatively  insensitive  to  J,  the  number  of  angular  grid  points,  for  the  chosen  values  of 
J .  The  number  of  time  steps  must  be  chosen  so  that  the  scheme  is  stable.  No  attempt 
was  made  to  determine  the  stability  condition  for  this  example.  Because  the  accuracy  of 
the  scheme  depends  on  the  three  parameters  /.  J,  and  K  it  is  difficult  to  discern  the  order 
of  accuracy  of  the  scheme. 


Tests  were  made  using  a  non-uniform  radial  grid  for  this  test  case.  We  found  that 
in  all  cases  it  degraded  the  accuracy.  Further  study  is  needed  to  determine  if  non-uniform 
grids  can  be  used  to  give  good  accuracy  and  less  restrictive  stability  limitations  on  the 
time  step. 

Table  2.  Wave  Equation  Results 


I  J  K 

ll«errll 

^err 

11  12  40 

.2050(-3) 

.351l(-3) 

21  12  80 

.5314(-4) 

.9447(-4) 

21  16  100 

.5294 (-4) 

.9525(-4) 

21  20  120 

.5397(-4) 

.9727(-4) 

41  12  160 

.1595(-4) 

.2998 (-4) 

41  16  180 

.1623(-4) 

.3117(-4) 

41  20  320 

.1662(-4) 

.3177(-4) 

The  second  test  problem  was  to  solve  the  elliptic  equation 

V^u  -  c(a:,y)u  =  0.  (7.5) 

on  the  unit  disk  with  u  specified  on  the  boundary.  The  exact  solution  was  given  by 

u(x,y)  =  exp((z  -  0.l)(y  -  0.5))  (7.6) 

with 

c(i,y)  =  (i  -  0.1)*  +  (y  -  0.5)*.  (7.7) 

The  equation  (7.1)  was  approximated  using  the  left-hand  side  of  (3.3)  for  the  lapla- 
cian  with  the  Fourier  method  being  used  to  approximate  the  derivatives  with  respect  to 
<i>.  The  solutions  were  obtained  with  the  LSOR  method  discussed  at  the  end  of  section  5. 
If  finite  difference  methods  are  used  in  the  angular  variable,  then  a  dirwt  solver  such  as 
that  of  Swarztrauber  and  Sweet  |1973)  can  be  used. 
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The  results  of  several  test  runs  are  displayed  in  Table  3.  The  number  of  grid  points 
is  given  along  with  the  iteration  parameter  u  and  the  tolerance  on  the  updates.  The 
iterative  procedure  was  stopped  when 

||u”+i  -  tt"||/u;  <  tol.  (7.8) 

The  number  of  iterations  required  for  convergence  is  seen  to  be  dependent  on  the  number 
of  radial  grid  points  and  not  on  the  number  of  angular  grid  points.  The  accuracy  at  the 
origin  is  seen  to  be  relatively  independent  of  the  value  of  J,  as  is  expected  from  the  analysis 
of  section  5.  This  was  also  noted  for  test  problem  1.  The  norm  of  the  error  is,  however, 
dependent  on  J  for  small  values  of  J.  If  J  is  sufficiently  large  then  the  second-order 
accuracy  of  the  scheme  is  seen. 

A  center  formula  similar  to  (3.4)  gave  results  comparable  to  those  obtained  by  (2.6). 
The  results  shown  were  obtained  by  using  center  formula  (2.6). 

TABLE  3.  Poisson  equation  results 


M 


r.VwV-' 


.y 


I  2  ^  1  W*  .. 

^  Ar*  rf  d<t>^  ^ 


(7.13) 


-  ^5 - ^-/.  ^75  j 


1  ^’■'.y  <1  ^  2  I 


rf  3<^* 


_  ILi  + 

rf  rf  3^  rd<f> 


=  u^  •  +  Au^ 

•li  *1^  ^ 


=  vr  +  Avr 


The  pressure  was  updated  by 


..+1  =  .  /l(^-n<^;^u-^-l<-^^) 

^  2Ar 

_  “i-Hj  ■*“»,/  ^  "t-2.j  ^^1,; 

6Ar  r, 


(7.14) 


where  7  is  an  iteration  parameter,  as  described  in  Strikwerda  j 1984b)  The  third-order 
differences  with  respect  to  7  in  (7.12)  and  (7.14)  are  necessary  to  preserve  the  regularity 
of  the  scheme  and  hence  the  smoothness  of  the  solution,  e.g.  Bube  and  Strikwerda  [1983] 
and  Strikwerda  [1984a].  The  derivatives  with  respect  to  ^  which  are  marked  with  a  carat 
in  (7.13)  and  (7.14)  are  computed  as  in  the  Fourier  method  with  the  addition  of  the  term 


(7.15) 


where  the  plus  sign  is  used  in  (7.14)  and  the  minus  sign  in  (7.13).These  terms  are  included 
to  ensure  the  regularity  of  the  scheme  and  hence  the  smoothness  of  the  solution.  Without 
terms  such  as  these  the  solution  would  contain  Fourier  modes  with  wavelength  2A^  of 
sufficient  amplitude  to  affect  the  accuracy  of  the  solution. 

The  results  of  test  problem  3  are  displayed  in  Tables  4  and  5.  In  Table  4  the 
norms  of  the  error  are  displayed  for  the  velocity  components  and  the  pressure.  That  is, 


7-1 J-l 


“erril  =  “  Z!  IZ  l«(’‘«7<^>)  "  Ui,;|*ri Ar A<> 


1=1 j=0 


where  the  initial  factor  of  is  included  to  normalize  by  the  area  of  the  disc.  The  error 
for  V  is  computed  similarly.  The  expression  u(ri,^y)  is  the  exact  solution  evaluated  at 
{ri,<l>j)  and  Uij  is  the  computed  solution  at  that  grid  point. 


Because  the  pressure  is  defined  only  to  within  an  additive  constant  we  use 


I  J-i 


Perrll  =  “  HZl  I  “  PiJ  “  ^P  I*  £.r,ArA(^ 


, 1=1  >=0 


(7.16) 


where 


if  t  =  0  or 
otherwise, 


/; 


as  required  by  the  trapezoid  rule  (6.7),  and  Ap  is  the  average  value  of  p(r,-,^,)  -  p,j 
computed  over  the  disc,  i.e. 


Ap  =  -  ^  ^(p(r„ <l>j)  -  p.j)e,rj Ar A(^.  (7.17) 

i=i  i=o 


Table  4.  Norm  Errors  for  Problem  3 


I  J 

Ikerrll 

l|t;.rrll 

llPerrll 

Sh 

11  16 

1.1(.3) 

9.1(.4) 

.34 

-1.5(-3) 

11  24 

3.9(-5) 

3.3(-5) 

3.2(-2) 

-5.9(-5) 

21  32 

2.l(-6) 

2.1(-6) 

4.0(-3) 

-2.3(-6) 

41  40 

l.l(-7) 

1.2(-7) 

3.6(-4) 

-1.5(.7) 

Also  displayed  in  Table  4  is  the  value  of  6^  which  is  the  average  of  the  discrete 
approximation  to  the  divergence  of  the  velocity  field.  The  finite  difference  and  Fourier 
scheme  does  not  enforce  the  condition 


V  •  u  =  0 

on  the  discrete  solution,  rather  the  iterative  method  converges  to  a  solution  with 

Vfc  •  u,.,  =  6„  (7.18) 


where  6^  is  the  average  of  the  left-hand  side  of  (7.18).  Thus  the  value  of  6h  is  an 


a  posteriori  indicator  of  the  accuracy  of  the  discrete  solution.  As  seen  in  Table  4  the 
numerical  method  gives  very  good  solutions  to  the  Stokes  equation. 

Table  5.  Errors  at  the  center  for  Problem  3 


Table  6  gives  the  iteration  parameters  and  the  resulting  number  of  iterations  for 
each  of  the  cases  reported  in  Tables  4  and  5.  The  iterative  method  was  considered  to  have 
converged  when  the  successive  changes  in  u  and  v  were  less  than  than  tol  times  u;  and 
when  the  changes  in  p  deviated  from  its  average  value  by  less  than  tol  times  'y.  That  is, 
when  the  value  of 


I  J-1 


1/2 


Vi  =  Oj=0  ) 

was  less  than  t(A  times  'y  the  solution  was  considered  converged. 


(7.19) 


The  values  of  the  expressions  (7.16)  and  (7.19)  was  computed  in  one  pass  through 
the  data  by  the  following  modification  of  West’s  algorithm  [West  (1979)].  To  compute  the 
value  of  t  j  I 

Q=  ''£(X„-Tl)'a„ 


■j=o 


where 


i.j-\  t.J-\ 

^  =  E  E  “o 

|,J=0  *.>=0 


initialize  with 


>1=0 
Q  =  0 
X  =  0 


Then  for  t  =  0  to  7,  and  j  =  0  to  J  -  1, 


Q 

X 

A 


Q  +_Aai,{Xii  -  Xy/{A  +  an) 
(AX  +  anXn)l{A  +  a.,) 


i4  +  a,j. 


(7.20) 


Thus,  to  compute  the  expression  (7.19)  we  have 


and 


a,j  = 


r,ArA<^,  i  ^  I 
|rjArA<^,  i  =  1,0. 


This  algorithm  makes  the  computation  of  expressions  (7.16)  and  (7.19)  only  slightly 
more  difficult  than  the  computation  of  the  usual  norm. 


Conclusions 


The  results  of  the  test  problems  show  that  the  methods  presented  in  this  paper  can 
be  used  to  compute  accurate  solutions  to  equations  on  domains  with  polar  grids.  The  basic 
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- v: 


formulas  can  be  used  with  most  numerical  procedures.  The  use  of  the  Fourier  method, 
while  not  essential  to  the  center  formulas,  is  very  convenient  and  efficient  to  use  with  polar 
grids. 
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